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ABSTRACT 

We carry out a series of high resolution (1024 x 1024) hydrodynamical simu- 
lations to investigate the orbital evolution of Jupiter and Saturn embedded in a 
gaseous protostellar disk. Our work extends the results in the classical papers of 
Masset & Snellgrove (2001) and Morbidelli & Crida (2007) by exploring various 
surface density profiles (a), where a cx r~ a . The stability of the mean motion res- 
onances(MMRs) caused by the convergent migration of the two planets is studied 
as well. Our results show that:(l) The gap formation process of Saturn is greatly 
delayed by the tidal perturbation of Jupiter. These perturbations cause inward 
or outward runaway migration of Saturn, depending on the density profiles on 
the disk. (2) The convergent migration rate increases as a increases and the type 
of MMRs depends on a as well. When < a < 1, the convergent migration 
speed of Jupiter and Saturn is relatively slow, thus they are trapped into 2:1 
MMR. When a > 4/3, Saturn passes through the 2 : 1 MMR with Jupiter and 
is captured into the 3 : 2 MMR. (3) The 3 : 2 MMR turns out to be unstable 
when the eccentricity of Saturn (e s ) increases too high. The critical value above 
which instability will set in is e s ~ 0.15. We also observe that the two planets are 
trapped into 2 : 1 MMR after the break of 3 : 2 MMR. This process may provide 
useful information for the formation of orbital configuration between Jupiter and 
Saturn in the Solar System. 

Subject headings: planetary systems:formation, planetary systems:protoplanetary 
disks, solar systems:formation 



1. Introduction 

Planet-planet interaction within an environment of gas disk is an important procedure 
that may account for the initial conditions of multiple planet system after the depletion of 
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gas disk, and thus affects their final orbital configurations. One of the notable configuration 
that two planets may achieve is the mean motion resonance (MMR). For example, the 
crossing of 2:1 MMR between Jupiter a nd Saturn was proposed to account for the later 
heavy bombardments of the Solar system (jTsiganis et al. 1120051 ; iGomes et al. 1120051 ) . MMRs 
are also common in the recently detected exoplanets. Among the ~ 45 detected multiple 
exoplanet systems, at least 7 planet pairs are believed to be trapped in low order MMRs 
(see Table [1] for a list). 

According to the general theory of disk-planet interaction, a single planet embedded 
in a gaseous disk may undergo various types of migration. For a planet with mass smaller 
than several Earth masses (M®), the angular momentum exchange between it and the gas 
disk will cause a net mo mentum loss on the planet and results in a fast orbit decay, which is 
called type I migration (lGoldreich fc Tremaine 1 11979c IWard 1119971 ). For a planet with mass 
comparable to Jupiter, it opens a gap around its orbit. Through the planet, the angular 
momentum exchange between the outer and inner part of the disk balances each other. This 
effect locks the planet and forces it to move as a part of the disk, at the viscous evolution 
timescale. This is called the type II migration (ILin fc Papaloizou II1986I ). Moderate mass 
planets(comparable to Saturn) may undergo very fast migration, which is believed to be 
caused by the large corotation torque risen from the perturbed coorbital zone of the planet. 
It is now named r unaway migration or type II I migration with timescale of several tens of 
orbit periods only ljMasset fc Papaloizou 1 120031 ) . 



Due to the different migration rates of two planets embedded together in a disk, MMRs 
may be established through the relative convergent migration. For example, when the outer 
planet is less massive than the inner one and undergoes type I or III migration, it will prob- 
ably catch up the more massive inner one which undergoes type II migration, even they are 
both migrating at the same direction. So, in multiple planet systems with heavier planet 
locating at the inner side, planets may be easily trapped into the MMRs. Many construc- 
tive studies have been made to investigate this scenario , eit her by three-body simu 



with prescribed disk effects, e.g., ISnellgrove et al. I (120011) and iNelson fc Papa 



oizou 



at ions 



(2002), 



or through h ydro dynamic simulations, e.g ., IPapaloizou fc Szuszkiewicz I ( 120051 ). Kiev et al. 
feooi l2005h and IPierens fc Nelson I pjol). 



Masset fc Papaloizou I ( 120031 ) investigated a case where the Jupiter(inner)-Saturn(outer) 



pair embedded in a protostellar disk. In that case, Saturn migrates inward very fast (type III 
migration) and is then captured into 3 : 2 MMR by Jupiter which undergoes slow type II mi- 
gration. As soon as the resonance is well established, the two planets reverse their migration 
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and move outward together, preserving their resonance. iMorbidelli fc Crida I (120071 ) con- 
firmed these results by using a more reliable code that describes the global viscous evolution 
of the disk. They considered a wide set of initial conditions and found the 3 : 2 MMR is a 
robust outcome of Jupiter-Saturn pair, which is compatible to the requirement of a compact 
initial configuration of Jupiter and Saturn (with period ratio slightly less than 2). After gas 
depletion, the subsequent 2:1 MMR crossing under th e interaction with planetesimal disk 
may achieve the present con figuration of Solar system tjTsiganis et al. Il2005t iGomes et al. 



20051 : IMorbidelli et al. 1120051 ). 



Morbidelli fc Crida I ( 120071 ) also showed that, the common migration rate of Jupiter and 



Saturn depends on the viscosity of the disk, as well as vertical scale height (H/r) of the gas 
disk. At H/r = 0.05 — 0.06, Jupiter and Sa turn seem to be in a quas i-stationary configuration 
without migration. With such conditions, IMorbidelli et al. I (120071 ) found that, the Jupiter- 
Saturn pair could maintain the 3 : 2 MMR for at least 1500 Jovian orbital periods when the 
gas disk dissipating slowly and smoothly. And we note that, through their evolutions, the 
eccentricities of Jupiter and Saturn stay at relatively low values(e < 0.05). 



In this paper, following the work of iMasset fc Snellgrove I ( 120011 ) and IMorbidelli fc Crida 



( 120071 ). we investigate the orbital evolution of Jupiter-Saturn pair embedded in gas disk with 
different slope of surface density (a), i.e., a = — In (a)/ In (r). We will show that, under suit- 
able disk condition, a > 4/3 (or < 1), Jupiter and Saturn will undergo outward (or inward) 
migration after trapped into 3 : 2 (or 2:1, respectively) MMRs. Thus the surface density 
profile plays an important role on determining the type of resonance and their consequent 
migrations. Also we will show that, under some circumstances, the eccentricities of Jupiter 
and Saturn will be excited up to 0.15 along the migration. The eccentricity excitation and 
subsequent stability of the MMRs between Jupiter and Saturn are discussed. Such a study 
help us to reveal the orbital architecture formation of our Solar system. 

Our paper is organized as follows: in Section 2 we describe our model including the 
numerical methods and the initial settings. In Section 3 we present our results as well as the 
analysis. The conclusions and discussions will be placed at Section 4. 



2. Model and Numerical Set up 

2.1. Physical model 

Following conventional procedures, we simulate the full dynamical interaction of a sys- 
tem includes a solar type star, two giant protoplanets and a 2-dimensional (2D) gas disk. 
The star is fixed at the origin of the system with both the planets and disk surrounding it. 
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We denote the inner and outer planets with subscript 1 and 2, respectively. 

For numerical convenience, the gravitational constant G is set to be 1. The Solar mass 
(M Q ) and the initial semi-major axis of inner planet (i?i = 5.2 AU) are set as the units of 
mass and length, respectively. Then the unit of time is j^Pw, where P±o the initial period 
of inner orbit. The mass of central star is set to be 1 Solar mass (M* = 1M Q ). 

We solve the continuity and momentum equations of the gas disk in 2-D space, neglecting 
the self-gravity of the gas. In order to properly describe the global viscous evolution of the 
whole disk and to avoid the annoying inner boundary, we employ a Cartesian grid. The 
vertically averaged continuity equation for the gas is given by 

da d(au x ) d(au y ) _ 

dt + dx + dy ' 1 ' 

where a is the surface density. The equations of motion in the Cartesian coordinates are 

d(au x ) | d(au 2 x ) | d(au x u y ) _ dP ^<9$ | dQ xx | dQ xy 
dt dx dy dx dx dx dy 

djauy) | d(au x u y ) | d(au 2 y ) _ dP ^ | dQ xy | dQ yy 
dt dx dy dy dy dx dy 

where P is the pressure and $ is the gravity potential of the whole system. <£> includes the 
softened potential of the central star softened potential of the giant planets (^1,^2) 
and the indirect potential $i n( j rises from the acceleration of the origin, which is caused by 
the planets and the gas disk. We adopt a softened gravity potential of the central star to 
avoid the singularity at the very center of the computational domain, 

®s = =. (4) 

a/x 2 + y 2 + el 

The potential of each planet is also softened, 

%, = - , GMp ' 1 (< = 1,2), (5) 

sj(x- x^i) 2 + (y- y P: ,i) 2 + e 2 p 

where e* and e p are the soften length to the central star and planets respectively. We test 
several values of e* and set it to be 0.05 for the balance of efficiency and accuracy. The e p is 
set to be O.QH, where H is the scale height of the gas disk. 
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The viscosity of gas is carried out by solving the stress tensors Q xx , Q yy and Q xy 
explicitly: 

where rj = av and 1/ is the dynamical viscous coefficient of the gas, which is assumed to be 
constant all over the disk. 

We assume the disk gas has a polytropic equation of state: 

P = Ka\ (9) 
where 7 = |. K is a constant that makes the following equation satisfied at r = 1, 



Cs(r)= vl^ = (7) ^p'^ 1 ' ( 10 ) 

where c s (r) is the speed of sound and v^ ep (r) is the Keplerian velocity of the gas. According to 
our settings and units, K = 0.1349. To focus on the effects of the surface density distribution, 
we fix the disk viscosity v = 10~ 6 and the disk aspect ratio H/r = 0.04. The planetary 
accretion and self-gravity of gas disk are neglected to reduce variables and to improve the 
computational efficiency. 



2.2. Mesh configuration and computational domain 

The Antares code we have developed is adopted in the simulations. It is a 2D Go- 
dunov code based on the exact Riemann solution for isothermal or polytropic gas, fea- 
tured wit h non-reflecting bo undary conditions. The detail of this method was described 



elsewhere (j Yuan fc Yen II2005I ). Since we adopt the Cartesian grids, the computational do- 
main is also a square. Figure [1] shows the computational domain. To get the proper gravity 
potential of a disk, we add an circular area outside the computational domain. Since it is 
far from the interested place, this circular area stays at the initial condition during the sim- 
ulations and its gravity potential is pre-calculated. The computational domain (gray area) 
is divided by a Cartesian mesh. The resolution of this mesh is a crucial issue for hydro- 
dynamics simulations, since poor resolution may introduce non-physical effects that affects 
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the migration of planets greatly. So we employ a high resolution: N x x N y = 1024 x 1024. 
The real computational time for a single 2000-orbits run is 3 — 4 weeks, in a full parallelized 
64-cpus cluster. This almost reaches the upper limit of normal computational ability, and 
the higher resolution may be achieved by the nested grids or other technics. 

The reasons a nd advantages of ou r choosing of Cartesian grids had been discussed in 
our previous paper ( jZhang et al. 1 120081 ) . For the orbits of two planets we integrate a three- 
body problem associate with the potential of gas disk by adopting a 8th-order Runge-Kutta 
integrator. The global time step is set as the minimum of the hydrodynamical and the orbital 
integration part. 



2.3. Corner-Transport Correction 

Although it is easy to implement, the Cartesian grids have a disadvantage to simulate 
a circularly orbiting gas disk, where the physical stream lines are neither parallel nor per- 
pendicular to either of the grid lines in x- and y-directions. The flux is evaluated at each 
interface at which the velocity of gas projects to the coordinate directions. When the grid 
resolution is low, the conservative law may break along the grid lines(see Figure [2]) and 
much dissipation arise. It is a kind of numerical viscosity which reaches maximum at the 
diagonal area of the computational domain. To make this numerical viscosity negligible, the 
resolution of the grids usually needs to be very high. Instead, we adopt the CTU(Corner- 
Transport Upwind) method to minimize this grid effect at the resolution that we can achieve 
so far. 

Following the CTU method, we add correcting terms to the flux at each interface of all 
cells. The value of this correcting flux depends on the angle between the direction of the 
physical velocity and the grid lines. As Figure [2] shows, by assuming u x > and u y > 0, 
the cell is affected by an additional flux comes from cell (i— l,j — 1). The cell averaged 
value Q™j, for example, is modified by a term, 

*^«B_ y -«.„-,), (ID 

where ^u x u v (At) 2 is the area of the small triangular portion moving into cell + 1) and 
Ax Ay is the area of the cell in which the jump Qfj — is averaged. The corresponding 
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fluxes evaluated at the four interfaces of cell now read 



F i^-l^v(Qi,3-Qi,3^ 



^+5 J 



(12) 



'« 2 



1 At 



U X Uy{Qij Qi-lj)- 



The details and im 



cussed by lColellal ( Il990h . 



2 A// " 

olementations of CTU method for conservation laws have been well dis- 



2.4. Comparison with other codes 

To ensure the reliability of our code, we performed a series of comparisons with other 
representative codes, e.g. FARGO. We examined the gap opening by Jupiter in a 2D 
gas disk. The phys i cal se tups and initial conditions are the same with that adopted by 



de Val-Borro et al. I ( 120061 ). except that our radial domain is from to 2.5 instead of [0.4, 2.5]. 
The grid resolution is 640 x 640 in the spatial range of [x min , x max } x [y xmin , y max \ = [-2.5, 2.5] x 
[—2.5,2.5]. Following their descriptions, we focus on the density contours of the gap, the 
evolution of density profiles, the evolution of total mass and the torques exerted on Jupiter. 
All the comparisons are performed in both inviscid and viscous disk, where the dynamical 
viscous coefficient is set to be v = and v = 10~ 5 respectively. The results with which we 
compared ourselves are obtained from the wetQ which maintained by de Val-Borro. 

Figure [3] shows the density contours after 100 orbits for the inviscid simulation. Two 
shocks are observed in our simulation (Antares): the primary one starts from the planet's lo- 
cation and the secondary one starts near the L5 point. We find the pattern of the spiral arms 
are similar with that of the other codes although the pitch angle of the primary arm(outside 
the orbit of planet) is a bit higher in our simulation, which is probably caused by the relative 
large sound speed their. Since the exact Riemann solution is not valid for locally isothermal 
gas, we adopt the full isothermal equation of state in this comparison tests instead and set 
the sound speed c s = — Vk ep \ r =i, where H is the height of disk and Vk ep is the Keplerian orbit 
speed. In the gap, there are two symmetric density enhancements locate close to the L 4 
and L5 points at azimuthal distance A<ft = ±7r/3 from the planet's location. This is in good 
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agreement with the theoretical prediction and the results of FARGO. We also notice there 
is a density bump, which indicates the congregation of vortensity, orbiting along the outer 
edge of the gap, whic h is observed in many oth er codes as well, e.g. the results of FARGO 



and NIRVANA- GDA (jde Val-Borro et al. 1 120061 ) . In the viscous simulation, the gap opened 



by the planet is narrower and smoother (Figure H]). The density enhancements seen at the 
Lagrangian points inside the gap in the inviscid simulation are dissipated, so is the density 
bump at the outer edge of the gap. Our simulation presents the same results with the other 
codes and shows more detailed structures within the gap. 

Figure [5] and [6] show the density profiles at different times in the inviscid and viscous 
simulations respectively. The width and depth of the gap in our simulations are in good agree- 
ment with those in FARGO. Our code also present the proper diffusion effects — decreasing 
on both the width and depth of the gap in viscous simulation. The major difference is the 
density on the inner disk: the surface density on the inner disk increases to higher value in 
our simulation than that in FARGO. This mainly comes from the different treatments at the 
center of disk. 

As the Cartesian grids are adopted, we need not introduce any inner boundary at 
the center of disk but just let the gas accumulate there. The competition between the 
increasing pressure, the dissipation of gas and the gravity of central star(there also could be 
the accretion of the central star which is not included in this comparison simulation) will 
result in an equilibrium and maintain an inner structure naturally (whose scale is around 
r < 0.1 and changes with time). In a real proto-stellar disk, there probably exists an inner 
boundary near the center of disk, however it should be maintained by some equilibriums, 
e.g. the evaporation and the refilling of gas, and should move inward or outward according 
to the changes of local situation, e.g. the enhancement of the density, instead of a fixed 
boundary. And further more, a full inner disk(includes the part r < 0.4) plays a great role 
on the dynamical interaction between the giant planet an d the global disk, and thus can not 



be ignored without careful treatments fjCrida et al. 1120081 ). 



Compared to the codes adopt absorbing or open inner boundary, our code maintains 
relative high level of the total mass in the disk during the simulation. After 200 orbits 
evolution, the total mass in our simulation decreases only ~ 2.5% for the inviscid case and 
~ 3.5% for the viscous case. While most results of the other codes are at the range of 
2% — 9% for the inviscid simulations and 5% — 8% for the viscous simulations (see Figure 

ED. 

The comparisons of torques are showed in Figure [8] (inviscid) and Figure M (viscous). 
Our results are coherent well with FARGO both qualitatively and quantitatively. The aver- 
aged total torque between 175 — 200 is —2.5 x 10 -5 in the inviscid simulation and —6.6 x 10~ 5 
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in the vis cous simulation. Bot h of th is two values are around the average of the results pre- 
sented by Ide Val-Borro et al. I (120061 ) . 

According to the above comparisons, we conclude that our code — Antares — is reliable 
for simulating planet-disc interaction in Cartesian grids. By adopting the CTU method 
and high resolution mesh, the numerical dissipation (or grid effects) is negligible in our 
simulations. 



2.5. Initial Conditions 



The surface density of the gas disk (a) varies as a function of its radius r. As shown 
in Table [2, we adopt several different initial distributions: ero, o"oe _r2//53 , oor -1 / 2 and so 
on. <T is set to be 0.0006 in our units, which corresponds to a height-integrated surface 
density ~ 200 g/cm 2 . The angular velocity of the gas ug = rQ g is slightly different from the 
Keplerian velocity since the flow is in a centrifugal balance with both the softened gravity of 
the star and the gas pressure which raises from the distribution of the surface density cr(r). 
The initial radial velocity of gas is set to be 0. These initial conditions of the disk do not 
take into account of the gravitational perturbation by the planets. 

To set up a dynamical equilibrium in which the orbit of a planet is circular and the 
stream lines are closed, we adopt an negligible initial mass for each of the planets (3 x 10~ 7 
or equivalently 0.1M e ). And at the very first 200 fixed circular orbits, both growth rates 
of the planets are specified to be ~ 3.5% per orbit until they achieve the mass of Jupiter 
and Saturn, respectively. With this "quiet-start" prescription, the planets gain their masses 
through adiabatic growth so that the disk has enough time to make a smooth response. The 
releasing consequence of the two planets is a complicate issue, since it directly relates to the 
formation consequence of a multiple planet system. We choose to release two planets at the 
same time, so the pre-formed planet should stay at an circular orbit and wait for the later 
one. Although this process may introduce some inconsistences, it's the most suitable initial 
state to investigate the migration of a planet pair. 



3. Results 



3.1. Effects of disk's surface density 



We consider a configuration in which S aturn locates outside the orb it of Jupiter, which 
is similar to our Solar system. As shown in iMasset fc Snell grove I (j200l[ ). this configuration 
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usually leads to the convergent migration and resonance trapping in a gas disk. In this 
paper we focus on the effects of the disk's surface density. The slope of the surface density 
determines the magnitudes of torques exerted on the planets, and thus affects the speed of 
the convergent migration as well as the type and stability of MMRs. We assume the surface 
density of disk is a function of its radii only: o = aor~ a and obtain a series of density 
distributions by varying the value of a. As summarized by Tabled most of the typical a 
had been considered: for example, a flat disk a = a , a moderate steep disk a = o"o r 1 and 
some extreme steep ones, e.g. a = <7 r~ 5 / 3 . 

We start with a flat disk, where the initial surface density is <jq = 0.0006 which cor- 
responds to a height-integrated surface density ~ 200 g/cm 2 . Panels (al-a4) in Figure [TU1 
show the evolutions of the semi-major axes, eccentricities and resonant angles of the two 
giant planets. When we release them at t = 200Pi , Jupiter had almost opened a gap and 
begins a slow type II migration after a short transition period. While the situation is quite 
different for Saturn: the tidal torque generated by Jupiter keeps pushing the gas into the 
coorbital zone of Saturn, thus greatly delays the gap opening process of Saturn. As a re- 
sult, Saturn migrates inward under the Lindblad and corotation torques, at a speed much 
faster than that of Jupiter. After the convergent migration, Saturn is trapped into the 2 : 1 
MMR with Jupiter. The eccentricities of both planets are excited as soon as the MMR is 
established. 

Following the flat disk, we run a moderate case in which the surface density of disk 



■2 



i s set to be a ~ e 53 . This density profile is derived from the analysis of iGuillot et al. 



(120061 ) for an evolving disk. It is in fact very close to the flat disk in our computational 
region(r e (0,4)). Embedding in this kind of disk, the migration of Saturn is a little more 
oscillatory than that in a flat disk, and we found the two planets stop and slightly reverse 
their migration to outward after the y had been locked in t o the 2 : 1 MMR. These are in 



good agreement with the results of iMorbidelli &: Crida I (120071 ) at the similar parameter 



settings(the same disk aspect ratio, viscosity and surface density profile), as shown in the 
Panels (bl-b4) of Figure [KB 

When the disk has a steep density profile, the situation becomes more complex. Accord- 
ing to our simulations, different values of a lead to different rates of convergent migration, 
thus results in the trap of different MMRs. The common migration speed and direction also 
depend on the slope of disk density. We will state them in details as follows. 

First, the rate of convergent migration before the trapping of MMR increases as the 
density slope a increases, as shown in Panel (b) of Figure [TTJ After the moment of release 
{t = 200Pio), the coorbital zone of Saturn is perturbed heavily by the density waves generated 
by Jupiter. The migration of Saturn is thus dominated by the corotation torque, which 
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scales with the gradi ent of the potential vorticity within the vicinity of the planet in the 
linear approximationf lGoldreich &: Tremaine 111979c IWard Ill99ll . Il992l ): 



Tc oc a 



d\og{cr/B) 
dlogr 



(13) 



-3/2 



in a nearly Keplerian disk, 



where B = k 2 / (4Q) is the second Oort constant with B 
and a = oo r_a is the surface density. So the direction and magnitude of the corotation 
torque, Tc oc (| — a)r~ a , depends on a. In fact, when the planet is forced to migrate fast 
(e.g., scattered by other planets) or its coorbital zone is perturbed heavily (in our cases, for 
example, the density waves generated by Jupiter), the density gradient in its coorbital zone 
is not monotone with r and usually very complex. So the corotation torque exerted on it 
should be evaluated by the real angular momentum exchanged within the coorbital zone. 
Masset fc Papaloizou I (120031 ) found a relation between the dramatic migration rate and the 
coorbital mass deficit 5m. In a Keplerian disk it reads: 



-aQ{M p 



5m)a = AT LR 



Tra 2 5m. 
( 

3x,o 



(14) 



where Q is the angular velocity of the planetary Keplerian motion, M p is the planet mass, 
is the differential Lindblad torque and x s is the half-width of planet coorbital zone. A 
runaway migration occurs when 5m becomes comparable to the mass of planet M p , because 
large 5m leads to fast migration which breaks the assumption that a stays constant and 
results in an instabilityf lMasset fc Papaloizou 1120031 ) . In our simulations, the steeper density 
profile between the orbits of Jupiter and Saturn results in heavier density waves that affect 
the coorbital zone of Saturn, and therefore, generates greater 5m. Panel (a) in Figure ITTI 
shows that the mass variation in the coorbital zone of Saturn increases at larger a. 

On one hand, Saturn migrates inward faster at large a. On the other hand, according to 
the viscous evolution of disk, Jupiter migrates outward when a > 1/2 (details are presented 
after Equation H7|) . So Saturn and Jupiter migrate to each other faster when the surface 
density becomes steeper, see Panel (b) in Figure [TTJ As a result, the time period that 
needed before the trapping of MMR between the two planets is longer in flatter disk. For 
example, it will take 1000-1500 P\q in a disk with a < 1/2 (Figure [TO]) , while it needs only 
300-500 Pio in a disk with a > 1 /2 (Figure 

Second, The two planets may be trapped into different MMRs when the surface density 
slope a varies. The phenomenon that different types of MMRs t hat can be trapped during 
conver gent migration depends on the migration speed is revealed in lPapaloizou fc Szuszkiewicz 
( 120051 ). In our cases, when the convergent migration is relatively slow, i.e., a < 1, Saturn is 
trapped by the (p + 1) : p — 2 : 1 MMR of Jupiter(Figure [T0lfT2l) . While Figure [13] shows 
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the results with a steeper disk where a = 3/2. The Saturn passes through the 2 : 1 MMR 
and is trapped by the 3 : 2 MMR of the Jupiter. Large p would be achieved by increasing 
the slope of disk profile (a), so that the speed of convergent migration is increased. However, 
we do not observe any resonances with p > 2 in our simulations — p = 2 is the highest value 
in our results. In fact, a transition from 3 : 2 MMR to 2 : 1 MMR is observed when the 
surface density of disk is very steep(a = 5/3), see Figure [TH In this simulation, Saturn 
first passes through the position of 2 : 1 MMR with Jupiter under fast inward migration and 
is then locked into the 3 : 2 MMR with Jupiter. When the resonance established, the two 
planets migrate outward together, and at the mean time the migration of Saturn becomes 
unstable as its eccentricity keeps growing. After several hundred orbits, when the eccentric- 
ity of Saturn grows to e > 0.15, the 3 : 2 MMR breaks and Saturn is scattered outward. As 
soon as the resonance breaks, eccentricities of both the two planets are damped effectively 
by the gas disk and th en Saturn is captured by the 2 : 1 MMR of Jupiter. This result is 



consistent with that of iPierens fc Nelson I ( 120081 ) . and we find the MMR is not so stable for 



high eccentricities of planets when they are embedded in a steep disk(a > 1). 

Third, the common migration speed and direction of the planet pair in MMR varies with 
a. The speed of common inward migration of the planet pair slows down as a increases. 
When a > 1, the two planets reverse their migration to outward, and the migration speed 
increases as a increases, see Figure [151 This phenomenon can be understood as follows. In 
the case of one planet, Jupiter undergo type II migration in the viscous disk. Each unit gas 
ring suffers a viscous torque 

rdQ 

T u = 2nr 2 ua^P, (15) 
dr 

where v is the viscosity, a = cror~ a is the surface density of gas disk, f2 ~ r~ 3 / 2 is the angular 
velocity of Kerplerian motion. Assuming a constant v across the disk, the viscous torque 
Y v = —STTvaor 1 /^ . By further assuming the disk remains Kelperian flow and integrating 
the motion equation in the azimuthal direciton, the angular momentum transportation per 
unit mass is governed by following formula: 

dr 2nra dr ' ^ 
where A denotes th e local injection rate of a ngular momentum per unit mass into the disk 



gas from the planet (ILin fc Papaloizou Ill986l ). The effect of A vanishes when the planet is 



treated as a part of the disk in type II migration. Then we can obtain the radial movement 
of the gas under the effect of viscous torque: 

1 8T u /dr 1 _ x 

- -3z/( a)r . (17) 



2ura d(r 2 Q)/dr v 2 
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Since Jupiter now follows the movement of the gas, this in fact results in an inward(or 
outward) migration of the planet at a < 1/2 ( a > 1/2, respectively). It can be seen from 
Panels (al) and (bl) in Figure [T2l the initial stage (t = 200 — 500Pio) of Jupiter evolution 
(before it meets Saturn). 

The presence of Saturn complicates the situation. Since Jupiter is massive than Saturn, 
Saturn is in fact forced to migrate with Jupiter when a > 4/3 (Figure [TBI and [T4"l) . In the 
case of 1/2 < a < 1, our simulation shows that the presence of Saturn perturbs the previous 
outward migration of Jupiter, and makes it undergo slight inward migration(see Figure [T2"j) . 
The reason is that when Jupiter encounters Saturn which is under fast inward migration, its 
outward migration is slowed down. Then, the eccentricity of Jupiter is excited by the 2 : 1 
MMR. As the orbit becomes eccentric, Jupiter cuts into the inner part of the gas disk and 
this causes a negative corotation torque. Figure [TBI shows the evolution of the torques exert 
on Jupiter and the mass variation within the coorbital zone of it. It clearly shows that, after 
t = 700Pio, the large mass variation results in a negative corotation torque that reverses the 
migration of Jupiter. The moment t = 700Pio, is just the moment that the eccentricity of 
Jupiter exceeds 0.17 which is the critical value required by Jupiter to cut into the inner gas 
disk (obtained by set ae > 2.5Rhui, with Rrui the Hill radius of Jupiter). 

In this section we show that, for Jupiter and Saturn embedded in a gaseous disk, their 
convergent migration rate, type of MMRs and subsequent common migration depend on 
the density slope a of the gas disk. The direct reason for these effects is the different 
torques(direction, value and type) exert on the planets, which vary with a. 



3.2. Torque analysis 

To understand the phenomena shown in previous section more deeply, we present here 
some torque analysis based on the linear estimation. After the release of planet at t = 200Pio, 
the migration of Jupiter couples with the response of gas disk. Since the coorbital zone of 
Jupiter is effectively cleared, the angular momentum transportation between the inner and 
outer disk are mainly done by the Lindblad torques. Due to the lack of Lindblad torques 
expression for a disk under the perturbation of Jupiter, we estimate the differ ential torque 



from linear estimation. The mth-order Lindblad torque can be expressed as tjWard 1 11997 



Papaloizou fc Terquem 1120061 ): 

,lr_ sign(fi p - tt)n 2 a(r) ^ 2 



m 3^^r+F(i+4e 2 ) 

where Q, fl p are the angular velocity at the specific location of Lindblad resonance and the 
planets, respectively. £ is a function that ensures the cutoff of the Lindblad torques at large 
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m. This is naturally satisfied when the density vanishes in the gap around the planet. \I/ is 
the force function which reads: 

* = + 2m2( "-% m . (19) 

By a high order interpolation and averaging along azimuthal direction, we obtain the nu- 
merical density a(r) and angular velocity Q(r) distribution in the disk at certain time. Then 
we determine the resonance positions for each m < 80 as well as the Lindblad torque rises 
there through Equations (fT8"j) and (fT9l) . 

In a flat disk, we find that the result follows the analytic prediction that the outer 
Lindblad torque is always stronger, see Figure [TTJ However, the existence of Saturn weakens 
the outer Lindblad torques exerted on Jupiter by pushing gas outward away. And the inner 
disk is strengthened by the steep density profile when a > 0. Figure [TBI shows the differential 
Lindblad torque exerted on Jupiter at the beginning of simulation (t — 0), the moment of 
release (t = 200Pio) an d the moment when the common gap has well formed (lOOOPio), i n 
a disk where a = 3/2. At the beginning of simulation, when the disk is still unperturbed, 
the numerical results are fit well with the analytic results except at 25 < m < 40. With the 
time passing by, both of the inner and outer torques decrease and the position of maximum 
moves toward smaller m, which corresponds to the gap formation process. At the release 
moment, when the Saturn has already formed, the outer torque decreases more than the 
inner one does, thus the net torque changes to positive. This positive torque drives Jupiter 
to migrate outward. When the gaps of the two planets overlapped, the inner and outer 
torques of Jupiter become comparable and are much smaller than the initial value by an 
order of 2 — 3. Then, Jupiter is in type II migration effectively. 

However, the migration of Saturn is dominated by the corotation torque. As we have 
mentioned before, the coorbital zone of Saturn is always perturbed by Jupiter when they 
are approaching to each other — even after the common gap has formed(Figure [i~9]) . The 
situation is more serious in a steeper disk, because the density profile may amplify the 
density waves generated by Jupiter. Figure [20] shows the torque evolutions of both Jupiter 
and Saturn embedded in a gas disk where o ~ r - 3 / 2 . Although the amplitude of corotation 
torque is comparable to that of Lindblad torque exerted on the Jupiter, it oscillates around 
zero and its average effect vanishes. The outer Lindblad torque is weakened by Saturn and 
the net differential torque is slightly above zero. 

At the mean time, the corotation torque exerted on Saturn overwhelms the Lindblad 
torques. At the moment of release when T = 200Pio, the corotation torque is negative 
since the density gradient is not perturbed much and still maintain negative around Saturn. 
Then, as the two planets approaching to each other, the density gradient within the vicinity 
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of Saturn changes to positive and the corotation torque evolves to positive as well. At this 
moment, Saturn lies close to the outer edge of the common gap, see Figure [2TJ When 3 : 2 
MMR is established, the orbit of Saturn becomes more and more eccentric. This makes 
Saturn cut through the outer edge of the gap periodically, and generates periodical torques. 
When the eccentricity of Saturn reaches 0.15, it is scattered outward by the increasing 
corotation torque and the resonance breaks(Figure [13] and [Tj 



As shown by Figure [22] the mass variation within the coorbital zone of Saturn and 
the corotation torques exerts on it are strongly correlated. In the plot, the short period 
corresponds to the orbit motion of Saturn, as Saturn cuts into the outer disk in every orbit. 
And the long period is the libration time of Saturn which reads: 



2na 
*iib = "B— 

tirn 



i on 

2 T ~dr 



(20) 



where R co is the half width of the coorbital zone and usually equals to 2.5 times of the Hill 
radius of Saturn. The libration period at the position of Saturn is ~ 32P 10 , as shown 
in Figure 122] from T = IOOOPiq to IIOOPiq. This is also the pe r iod that the coorbital 



mass exchanges its angular momentum with the planet ([Ward Ill99ll ; iMasset fc Papaloizou 



20031 ). The mass variations are normalized by the unperturbed coorbital mass of each planet 
respectively. For Saturn, the peak value immediately before the onset of instability is about 
5 times greater than its unperturbed state, which is roughly 10 -4 in our unit, or 30% of the 
mass of Saturn. The angular momentum exchange between this part of coorbital mass and 
Saturn results in the corotation torque that dominates Saturn's migration. 

In this section we show that the different migration rates and directions of the two giant 
planets are the results of the combination of Lindblad and corotation torques, which depend 
on the density slope of the disk (a). The corotation torque exerted on Saturn weakens the 
stability of Saturn's orbit. The unstable migration is more serious in the disk with large a, 
see Figure [T2lfT^l And furthermore, the breaks of MMR are observed when a > 3/2 (Figure 
32] and [HJ). To analysis these breaks is helpful to understand the orbital evolution following 
the convergent migration of Jupiter and Saturn. The details and analysis are presented in 
the next section. 



3.3. Stability of MMRs 

MMRs are common results of the convergent migration of Jupiter and Saturn. In this 
section, we investigate the stability of them. As showed in Table [2] Saturn may be trapped 
by the 2 : 1 MMR of Jupiter when the disk is nearly flat and may reach the 3 : 2 MMR if 
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the surface density profile is steep, e.g. a > 1. Furthermore, the 3 : 2 MMR of the Jupiter 
and Saturn is not stable as their eccentricities keep growing, while the 2 : 1 MMR seems to 
be more robust even at relatively high eccentricities ~ 0.2(e.g., Figures [TU1 and [T2"j) . 

The instability of planets in the 2 : 1 or 3 : 2 MMRs could be caused by the overlap of two 
nearby resonances when the eccentricity is large enough. In the case that both two planets 
are with equal masses (m) and move in circular orbits, the overlap of ne arby resonances 



occurs when th e difference of their semi- major axes is smaller than a limit ( jWisdom 111980 
Gladmanlll993h . 

Aa 2 . m , 9/7 , . 

a dp 

where the relationship of n\jni = (p + l)/p = [(a + Aa)/a] 3 / 2 ~ 1 + 3/2(Aa/a) is used, 
with ni, «2 the mean motion of the two planets. Assuming m/M* = (0.27 — 0.95) x 10~ 3 for 
Jupiter and Saturn, the application of Equation (l2~Tj) give p > p min ~ 3 — 4 for the overlap 
between (p + 1) : p and (p + 2) : (p + 1) MMRs. 

However, when p < p m in their MMRs can also overlap if they are in high eccentric 
orbits. In the framework of restricted three-body problem, a planetesi mal embedd e d in a 



gas disk will be trapped into the outer resonance of a massive body. iKary et al. I (119931 ) 
estimate of the minimum separation between MMRs below which the planetesimals may 
undergo chaotic instability due to the MMR overlap: 

Aa 8nm p 2 3v di{ 1/9 

a ~ U 3M» 1 { 2v kcp )l ' 1 1 

where m p refers to the mass of planet and e is the orbital eccentricity of planetesimal, v dif is 
the difference between the gas velocity and the local circular Keplerian velocity v kep . This 
velocity difference is a result of the pressure rose by the slope of the surface density a in the 
disk, which makes the gas circle the central star at a sub-Keplerian velocity, and is referred 
as the drag effect of the gas disk. In a nearly flat disk, Vdit/v^p is usually tiny, for example, 
^difAVp = 0.002 at r = 1 when cr ~ e _r2//53 , while it increases to ~ 0.01 when a > 3/2. 
Use the approximation t h at 3/ 2(Aq/a) 1/p ~ l/(p+ 1) and the relation obtained by 



Weidenschilling fc Davis I (119851 ): 



^ ( ^p )1/2) (23) 

we can get the minimum eccentricity, above which the (p+ 1) : p and (p + 2) : (p+ 1) MMRs 
will cross, from Equation ( 122]) : 

3M,,2 1 „ 



Sirm p 3 p + 1 
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This gives e m in ~ 0.065 for p + 1 : p = 3 : 2 MMR, and ~ 0.49 for 2:1 for a Jupiter mass 
disturber. Although Equation (|24p is obtained from the restrictive three-body problem, 
our simulation show that it also applies to Jupiter-Saturn case. In fact, we find the 3 : 2 
MMR breaks as soon as the eccentricity of Saturn reaches 0.15 and larger a doesn't change 
this value but only makes instability happen earlier, see Figure [13] and [HJ The 2 : 1 MMR 
is much more robust, because it only requires e < 0.5 to maintain stable and the eccentricity 
excited by 2 : 1 always stay at e < 0.2, see Figures [TUl and [i~2l 

The other reason for the instability of MMR is the runaway migration induced by the 
corotation torque. We had shown that, a moderate planet will undergo runaway migra- 
tion(Type III migration) when its c oorbital zone is perturbed heavily by the density waves 



generated by the other giant planet ( IZhang et al. M2008I : IZhang fe Zhou 1120081 ) . It is believed 



that two giant planets would be trapped and migrate together when their gaps overlapped 
well. It's true that when the density distribution of disk is nearly flat, the convergent mi- 
gration is relatively slow and Saturn has enough time to clean its coorbital zone before it 
interacts with Jupiter directly. However, when the disk is steep, say a ~ r ~ 3 / 2 ? the conver- 
gent migration is so fast that the Jupiter would catch Saturn to the 3 : 2 MMR before the 
later one opens a clean gap, see Figure |23J At this time, Saturn is forced to migrate with 
Jupiter while the corotation torque still dominates its migration. As we can see from Figure 
[TBI the migration curve becomes very oscillatory as soon as Saturn is trapped by the MMR 
with Jupiter. 

In fact, eccentricity plays an essential role on triggering runaway migration even if a 
clean gap formed around the planet orbit. As we mentioned in previous section, as soon as 
the giant planet undergoes eccentric motion, it cuts through the edge of the gap periodically 
when ae > -R gap , where R gap is the width of the gap. Then, a periodic torque rises as a 
result of the replenishment of gas into the vicinity of the planet. This torque may change 
the direction of planet's migration (see Figure [TB"]) or kick the planet away quickly. A good 
approximation for the gap width is i? gap ~ 2.5i?Hiii, the same with the radius of the coorbital 
zone of a giant planet by assuming the viscosity is low, where Rmn is the Hill radius of this 
planet. In the case of Saturn, it requires e ~ 0.15 to allow Saturn cut into the disk. However, 
as Saturn lies outside the orbit of Jupiter, it can only touch the outer edge of the common 
gap, where the surface density gradient is positive and the corotation torque exerted on it is 
positive as well, see Figure [21] Thus, Saturn tends to be scattered outward at this kind of 
configuration. 

According to our results and analysis, we find that the instability of MMR is mainly 
due to the large eccentricity excited by the MMR itself. High eccentricity either leads to 
the overlap of MMRs or makes the planet cut into the disk and generates strong corotation 
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torque. As the ecce ntricity evolutions of the two planets locked in resonance vary with 
the type of resonance (jMichtchenko et al. 1120061 ). density slope (a) will affect the final orbit 
configuration of the planet pair. 



4. Conclusions and discussions 

A series of high resolution hydrodynamic simulations have been performed to investigate 
the orbital evolution of Jupiter and Saturn embedded in a protostellar disk. We focus on 
the effects of different surface density profiles of the gas disk where a oc r~ a . The instability 
of mean motion resonance caused by the convergent migration is also studied. According to 
the results and analysis of our simulation, we summarize our conclusions as follows, as well 
as the discussions and implications for these results. 

(1) The existence of Jupiter (massive planet) delays the gap formation process of Sat- 
urn(light planet) and generates great perturbations within the coorbital zone of Saturn. 
These perturbations result in the inward or outward runaway migration(type III migration) 
of Saturn, depending on a, the density slope of the gas disk. 

The effects of the pre-formed giants are very important for understanding the orbital 
evolution of a multiple planet system. To investigate these effects, we fix the very first 200Pio 
orbits of both the two planets when they are growing from 0.1 Earth mass to 1 Jupiter and 
1 Saturn mass respectively, then we release them at the same time. This, in fact, assumes 
that the Jupiter had formed before the Saturn did since the gap of Jupiter had well formed 
while the Saturn's hadn't yet. 

The orbital evolutions of the planets form later are greatly affected by the pre-formed 
giant ones. The first generation of giant planets stro ngly modified their n eighborhood by 



opening gaps at se veral Hill radius fr o m themselves (IBrvden et al. 1 120001 ) . as well as the 
planetesimal gaps Jzhou fc Linll2007h . IPierens fc Nelson I J2008T had shown that the light 



planets will be trapped at the edge of the giant planet's gap. This halt is because of the 
balance between the Lindblad torque and the corotation torque which rises at the density 
jump. This happens when the mass of the outer planet is low: ni2 < 20M®. For massive 
outer planets, i.e., m 2 ~ Mj, the strong tidal effect guarantees the gap formation process and 
slow type II migration could be expected. The situation becomes complex when the outer 
planet has a moderate mass, m 2 ~ Mg which is critical to open a gap. Its orbital evolution 
thus strongly depends on the initial conditions, for example, the initial density profile of the 
disk a. 



Zhang et al 



(120081 ) had shown that the existence of the gaseous disk expands the 
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interaction region of the proto-planets. As we have shown in the previous section, the tidal 
effect of Jupiter keeps pushing the gas away from its coorbital zone and replenish the gap of 
Saturn (see Figure [T91 and |2~TT) . This effect increases as a increases(Figure [TT]) and drives 
the orbital evolution of the system unstable. If the initial separation is large or the disk is 
nearly flat, clear gaps would form around the planets' orbit and gentle convergent migration 
could be expected. This will result in a less compact and more stable system. 

(2) The convergent migration rate is proportional to the surface density slope a. And, 
as a result, the types of MMRs depend on a as well. The two planets approach to each other 
gently when a < 1, and are locked into 2 : 1 MMR. When a > 1, the convergent migration 
is fast and the 3 : 2 MMR is reached. 

The convergent migration rate is one of the most critical issues to determine the con- 
sequential migration of the planet pair. Many facto rs may affect this rate, for example, the 



masses of the two planets (IPierens fc Nelson II2008I ). the viscosity and the disk aspect ratio 



( jMorbidelli fc Crida 1120071 ) . However, the essential factor which results in differential migra- 
tion is the various torques by which the planets are driven. These torques intensively relate 
to the surface density of the disk in which the planets embedded. For a low mass planet, 
the different ial Lindblad torque is negative and the value is proportional to the density slope 



a(a ~ r~ a ) (lWard Ill99lf ). The migration of massive planet is determined by the global dis- 
tribution of the angular momentum(relates to a) and the viscosity of the gaseous disk. The 
corotation torq ue relates to the vortensity g radient (relates to a as well) within the coorbital 



zone of planet ( Masset &: Papaloizou 20031 ). The analysis of torques by a semi-analytical 



method have been shown in the previous section, see Figure [ToTfTHl and 

To get the different convergent rate naturally, we choose a series of density profiles of 
the disk where a ~ r~ a . Most of the typical density profiles have been adopted, see Table 
[2J Our results show that the convergent migration rate increases as a increases, see Figure 
rm This is reasonable since the migration of Saturn is accelerated by the steep density 
slope while the migration of Jupiter turns outward when a > 1/2. The type of resonance 
is determined by how close the two planets could approach. When the disk is nearly flat, 
where a < 1, the 2 : 1 MMR is a robust outcome. While the 3 : 2 MMR is more favorite in 
the steep disk where a > 1. The disk with a density profile around a ~ 1 shows a transition 
phenomena and the high eccentricity of Jupiter reverses its outward migration to inward, 
see Figure [121 and I 



(3) The 3 : 2 MMR of the two planets is unstable when the eccentricity of Saturn 
becomes large enough in a steep disk where a > 3/2. We estimate that the critical value is 
e s ~ 0.15 with our settings. 
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The 3 : 2 MMR of Jupiter and Saturn is thought to be robust for many settings, but 
we find that this configuration may breaks down when the eccentricity of Saturn grows too 
high in a gaseous disk where the surface density gradient is pretty steep a > 3/2. In fact, 
the eccentricities will be exited as soon as the resonance established no matter the disk is 
nearly flat or very steep (Figures [TU1 and fT2lTH|) . However, in a steep disk, the situations are 
different. First, the convergent migration is much faster. This enables the two giant planets 
get much closer and as a result, the dynamical instability becomes possible, for example, 
the overlap of resonances. Second, caused by the fast convergent migration as well, the 
establishment of resonance occurs before a clear common gap formed. This will result in a 
strong corotation torque which may drive an unstable migration of Saturn, see Figure [13j 
Third, the steep gradient of density in fact amplifies the density waves propagating from 
Jupiter to Saturn and results in relatively heavier perturbation in coorbital zone of Saturn. 

In our simulations, the planet pair migrates outward when a > 1. The separation 
between them increases as they preserve the 3 : 2 MMR. Thus, Saturn is pushed outward 
further and further. At the mean time, the growing eccentricity makes Saturn cut into the 
outer disk deeper and deeper. The density gradient at the out edge of the common gap is 
positive and generates positive corotation torque that pushes Saturn outward further. When 
the eccentricity is high enough, Saturn is scattered outward. Our results show that the 
critical value is e > 0.15 ~ 0.2 for 3 : 2 MMR of Jupiter and Saturn. 

We also find that the onset of instability would be suspend when a decreases. In fact, 
in the case where a = 4/3, the two planets maintain 3 : 2 MMR over 2000Pio, see Figure 
[T5l So the long time stability could also be expected by choosing a proper a, and this still 
needs further simulations. The 2 : 1 MMR seems to be more stable for high eccentricities 
and we find the two planets could be re- locked into 2 : 1 MMR just after the break of 3 : 2 
MMR when they are both migrating outward, see Figure [14] 

If the initial density profile of our Solar nebular is relatively steep, e.g. a > 4/3, then 
the formation of the main configuration of our Solar system could be this way: (a) Jupiter 
and Saturn first migrate inward when their masses are still low. (b) By accreting gas from 
the nearby disk, they grow up quickly. While Jupiter grows much faster than Saturn does 
through the runaway accretion, and the mass difference results in a convergent migration 
between them, (c) Then the 3 : 2 MMR of Jupiter and Saturn is established and their 
migration turns outward with the resonance preserved, (d) When the eccentricity of Saturn 
becomes too high, the resonance breaks down and Saturn is scattered outward to the place 
near its present location or re-captured by the 2 : 1 MMR of Jupiter. Neptune and Uranus 
are also scattered out away at the mean time, (e) The eccentricities of both the inner rocky 
planets and the outer gas giant are damped effectively by the gas disk after the instability, (f) 
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After the dissipation of gas, the planets evolve to their present locations by the interaction 
with the planetesimal disk. Of course, many details are need to be addressed especially 
for the orbital ev o lution s of the inner rocky planets. However, compared to the results of 



Morbidelli et al. I ( 120071 ). our results suggest that the instability would happens before the 
dissipation of gas. The remaining gas may corresponds to the low eccentricities of the main 
planets in our Solar system. And since Jupiter migrates inward before its outward migration, 
the final location is not far from its birth place. 

We also notice that the excitation of Jupiter's eccentricity is previous than that of 
Saturn in 2 : 1 MMR and is laggard in 3 : 2 MMR. Since the eccentricity is a critical issue 
to guarantee the stability of Jupiter-Saturn pair, its evolution and constraints need to be 
addressed in details by considering the effect of the interacting disk. The results of a reverse 
orbital configuration of Saturn and Jupiter is in preparing as well. 
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Table 1. A short list of extra-solar planets in resonance. 



System 


1NO. 


f [dayj 


A/T cin ^ r A/7" 1 

ivi sini 


n \ ATT} 

a [AU\ 


e 


rteis. 


GJ 876 


C 


30.57 


0.56 


0.13 


0.2 


1 


(2:1) 


b 


60.13 


1.89 


0.21 


0.04 




HD 128311 


b 


458 


2.3 


1.08 


0.23 


2 


(2:1) 


c 


918 


3.1 


1.71 


0.22 




HD 73526 


b 


188 


2.90 


0.66 


0.19 


3 


(2:1) 


c 


377 


2.50 


1.05 


0.14 




HD 82943 


b 


217.9 


1.4 


0.74 


0.46 


4 


(2:1) 


c 


456.6 


1.78 


1.19 


0.36 




HD 160691 


d 


310.55 


0.52 


0.92 


0.067 


5 


(2:1) 


b 


643.25 


1.68 


1.5 


0.13 




HD 45364 


b 


227 


0.19 


0.15 


0.17 


6 


(3:2) 


c 


343 


0.66 


0.68 


0.09 




HD 60532 


b 


201 


3.15 


0.76 


0.27 


7 


(3:1) 


c 


605 


7.46 


1.58 


0.038 





Note. — References. (1) Marcy et al. 2001; (2) Vogt et al. 2005; (3) Tinney et al. 2006; 
(4) Lee et al. 2006; (5) Gozdziewski et al. 2007; (6) Correia et al. 2009; (7) Laskar & Correia 
2009. 

Table 2. A summary of our simulations. 



Case Configuration Direction a Resonance Instability 



1 


J-S 


inward 




2 ; 


1 


no 


2 


J-S 


inward 


a ~ e .w 


2 : 


: 1 


no 


3 


J-S 


inward 


i 

a ~ r 2 


2 : 


: 1 


no 


4 


J-S 


inward 


•2 

a ~ r 3 


2 : 


: 1 


no 


5 


J-S 


inward 




2 : 


: 1 


no 


6 


J-S 


outward 


4 

o ~ r 3 


3 : 


: 2 


no 


7 


J-S 


outward 


a 

a ~ r 2 


3 : 


: 2 


yes 


8 


J-S 


outward 


5 

a ~ r 3 


3:2^2:1 


yes 
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Fig. 1. — The computational domain is from -4 to 4 in x direction and from -4 to 4 in 
y direction (gray square). Surrounding it are four non-reflecting boundaries. Area outside 
the square is assumed to stay constant. We take the gravity comes from the whole round 
area(.R < 6) as a background potential which is a function of radius. 
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Fig. 2. — The real flux flows at an angle to the grid lines(assume u x > 0,u y > 0). (a) 
The new value Q^j 1 in cell should also be affected by the old value Q^j^ in cell 
(i — 1, j — 1). (b) The flux at four interfaces of cell need additional corrections — see 
the four dark triangles. 
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Fig. 3.— Density maps in logarithmic scale after 100 orbits for the inviscid simula- 
tions. Note that the computational domain of Antares is (r m i n ,r max ) = (0,2.5), instead 
of (0.4,2.5). The density range is —1.7 < log(o"/o"o) < 1. Antares adopts isother- 
mal equation of state(EOS) instead of the locally isothermal EOS, so the pitch angle 
is bit large than that in FARGO. The results of FARGO are obtained from the web: 
http://www.astro.su.se/groups/planets/comparison, which is maintained by de Val-Borro. 
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Fig. 4. — Density maps in logarithmic scale after 100 orbits for the viscous simulations. The 
viscous coefficient v = 10~ 5 . The density range is also —1.7 < log(cr/cr ) < 1. The results 
of FARGO are obtained from the web: |http : / / www. astro .su.se/ groups / planets / compar ison , 
which is maintained by de Val-Borro. 
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Jupiter Inviscid. 
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Fig. 5. — The normalized surface density profiles averaged azimuthally over 2tc after 100 
orbits for the inviscid simulations. Note that the computational domain of Antares is 
( r min, r max) = (0,2.5), instead of (0.4,2.5). The surface density at inner disk increases 
and keeps a high level when the inner open boundary is absent. The results of FARGO 



are obtained from the web: http: / /www.astro.su.se/group s/planets/comparison[ which is 
maintained by de Val-Borro. 
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Fig. 6.— The normalized surface density profiles averaged azimuthally over 2tt af- 
ter 100 orbits for the viscous simulations. The viscous coefficient v = 10 -5 . Antares 
shows the proper behavior of gas under the increasing dissipation effect — the gap be- 
comes narrower and shallower. The results of FARGO are obtained from the web: 



http://www.astro.su.se/groups/planets/comparison which is maintained by de Val-Borro. 
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Fig. 7. — The evolution of total mass contained in the disk, (a) Total mass of disk calculated 
by Antares in inviscid and viscous simulations, (b) Total mass calculated by other codes in 
inviscid simulations. The legend details were presented in de Val-Borro (2006). (c) Total 
mass calculated by other codes in viscous simulations. See de Val-Borro (2006) again for 
legend details. 
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Fig. 8.— Comparison of the time-averaged torques exerted on Jupiter in 
inviscid simulations. The results of FARGO are obtained from the web: 
http : / / www. astro .su. se / groups / planets / comparison \ which is maintained by de Val-Borro. 



-32 - 




Fig. 9. — Comparison of the time- averaged torques exerted on Jupiter in viscous simulations. 
The viscosity coefficient is v — 10 -5 . The results of FARGO are obtained from the web: 



http://www.astro.su.se/groups/planets/comparison, which is maintained by de Val-Borro. 
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Fig. 10. — Orbital evolutions of Jupiter and Saturn embedded in a flat disk where a = a^r 
(Panels al-a4) and a nearly flat disk where a = cr e_r2 ^ 53 (Panels bl-b4), respectively, (al) 
and (bl) Evolutions of the semi-major axes of the planets. The two planets approach to each 
other gently in a flat disk or a nearly flat disk. The common migration stops or even reverses 
after the establishment of resonance. (a2) and (b2) Evolutions of the eccentricities of planets. 
(a3) and (a4) Evolution of the resonance angles:^! = 2X 2 — X± — w\ and # 2 = 2A 2 — Ai — zu 2 
of 2 : 1 MMR. (b3) and (b4) Evolution of the resonance angles: 9\ = 2X 2 — Ax — vj\ and 
6 2 = 2X 2 — Ai — -072 of 2 : 1 MMR. (A color version of this figure is available in the online 
journal.) 
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Fig. 11. — (a) Mass variation in the coorbital zone of Saturn during the convergent migration 
stage. It decreases when a decreases (a oc r~ a ) — from the top to the bottom in this panel, 
respectively. The radii of coorbital region is set to be 1.5Rhm and the mass variation 
is normalized by the initial mass contains in Saturn's coorbital region, (b) The different 
convergent migration rates with various a: it increases at large a. The convergent migration 
is halted when the two planets are locked into MMRs(2 : 1 for a < 1 and 3 : 2 for a > 1). 
Note that mass variation decreases as the common gap forms. In Panel (a), the large 
variations in the curve of a oc r -5 / 3 when T > 400Pio are the results of the increasing 
eccentricity of Saturn. This corresponds to the oscillations of the convergent migration(a = 
5/3) in Panel (b). (A color version of this figure is available in the online journal.) 
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Fig. 12. — Orbital evolutions of Jupiter and Saturn embedded in a slightly steep disk where 
a = ao r (Panels al-a4) and a steeper disk where a = uor^ 1 (Panels bl-b4). (al) and 
(bl) Evolution of the semi-major axes of the two planets. Jupiter first migrates outward 
when a > 1/2 (a oc r~ a ) and then reverses its migration to inward after locking into 2:1 
MMR with Saturn. (a2) and (b2) Evolutions of the eccentricities of the planets. Note that 
the Jupiter's eccentricity is excited heavily by the resonance. (a3) and (a4) Evolution of 
the resonance angles: Q\ = 2A 2 — X± — W\ and 82 = 2X2 — X± — W2 of 2 : 1 MMR. (b3) 
and (b4) Evolution of the resonance angles: Q\ = 2A 2 — Ai — w\ and 6 2 = 2A 2 — X\ — tu 2 
of 2 : 1 MMR. Due to the effects of the dissipative disk, the libration center changed from 
(61,82) = (0°,0°) asymmetric libration. (A color version of this figure is available in the 
online journal.) 
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Fig. 13. — Orbital evolutions of Jupiter and Saturn embedded in a disk where a = a r^ 3 ^ 2 . 
Saturn passes through the position of 2 : 1 MMR of Jupiter and then is catched by the 3 : 2 
MMR. Scattering happens at T = 1500Pio when e s > 0.15. (a) Evolutions of the semi- 
major axes of the two planets. Migration of Saturn becomes unstable after it is trapped into 
resonance with Jupiter, (b) Evolutions of the eccentricities of the planets. The excitation 
of Jupiter's eccentricity by the 3 : 2 MMR is behind that of Saturn, (c) Evolutions of the 
period ratio Ps/Pj = (0-2/ a i) 3 ^ 2 - (d) Evolutions of the resonance angle 6\ 2 — 2Ai — w 2 of 
3 : 2 MMR. The resonance becomes unstable as the eccentricities keep growing. (A color 
version of this figure is available in the online journal.) 
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Fig. 14. — Orbital evolutions of Jupiter and Saturn embedded in a disk where a = a r~ 5 / 3 . 
Saturn passes through the position of 2 : 1 MMR of Jupiter and then is catched by the 3 : 2 
MMR. Scattering happens at T = 8OOP10 when e s > 0.15. Eccentricities decrease rapidly 
after the break of resonance and then the Saturn is trapped by the 2 : 1 MMR of the Jupiter 
at T = 1400Pio • ( a ) Evolutions of the semi-major axes of the two planets, (b) Evolutions of 
the eccentricities of the planets. It clearly shows that the excitation of Jupiter's eccentricity 
is previous than that of Saturn in 2 : 1 MMR and is laggard in 3 : 2 MMR. (c) Evolutions 
of the period ratio Ps/Pj = (0-2/ a i) 3 ^ 2 - One can see the re-capture of 2 : 1 MMR. (dl) 
and (d2) Evolutions of the resonance angle: 9 = 3A 2 — 2Ai — ro 2 of 3 : 2 MMR(dl) and 
9 = 2X2 — X\ — w\ of 2 : 1 MMR(d2). (A color version of this figure is available in the online 
journal.) 
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Fig. 15. — Migrations of Jupiter and Saturn embedded in different disk where the density 
slope a varies from to 3/2. The two planets migrate inward when a < ~, and migrate 
outward when a > 1. It shows a transitional state when ^ < a < 1. (A color version of this 
figure is available in the online journal.) 
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Fig. 16. — (a) Evolutions of the semi-major axes of Jupiter and Saturn, (b) Evolutions 
of the inner, outer and net Lindblad torques as well as the corotation torque exerted on 
Jupiter. The net Lindblad torque is positive at the release moment (t = 200Pio), then a 
negative corotation torque rises and dominates Jupiter's migration(t > 700Pio). (c) Mass 
variation within the coorbital zone of Jupiter. The large oscillations after t = 700Pio result 
in the negative corotation torque which reverse the migration of Jupiter. The surface density 
profile is a = a^r in this simulation. (A color version of this figure is available in the online 
journal.) 
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Fig. 17. — The m-th inner and outer Lindblad torques exerted on Jupiter at different times, 
where |m| = 2 — 80. The density profile of the disk in which the two planets embedded is 
a = (Tor . From the top to the bottom, the time points are the initial moment, the release 
moment and the moment after common gap has formed, respectively. At initial moment, 
the torques are calculated as if Jupiter was embedded in an unperturbed disk. The solid 
lines denote the analytic results while the diamonds and circles denote the semi-analytic 
results. Note that the outer Lindblad torque is always larger than the inner one and the 
position of maximum value moves left as the gap becomes deeper and wider. The kinks in 
the numerical results indicate the sudden changes of density or angular velocity distribution 
of gas, or the failures of locating the exact positions of resonances( extrapolation value is 
adopted instead), which is mainly due to the lack of radial resolution of the numerical data. 
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Fig. 18. — The m-th(|m| = 2 — 80) inner and outer Lindblad torques exerted on Jupiter 
which embedded in a disk where a = o§r~^l 2 . From the top to the bottom, the time points 
are the initial moment, the release moment and the moment after common gap has formed, 
respectively. At initial moment, the torques are calculated as if Jupiter was embedded in an 
unperturbed disk. The solid lines denote the analytic results while the diamonds and circles 
denote the semi-analytic results. Note that the inner Lindblad torque becomes larger than 
the outer one at the release moment. Compared to Figure [171 one may find that the inner 
Lindblad torque almost remains the same level while the outer one decreases a lot. This 
could be explained by the mass reduction in the outer disk caused by Saturn. 
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Fig. 19. — This figure show the cross section of surface density at different moments. The 
initial density profile is a = o"o r_1 • One can see the overlapping process of two gaps. Note 
that the gap formation of Saturn is much delayed(T w 1000Pio)and the variation of density 
is acute at Saturn's vicinity. The coordinates of planets in the figure denote their positions 
(x axis) and the average gas density around them (y axis). (A color version of this figure is 
available in the online journal.) 
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Fig. 20. — (a) Evolution of the semi-major axes of Jupiter and Saturn, (b) Evolution of 
the inner/outer Lindblad torques and the corotation torque exerted on Saturn. It shows 
clearly that the migration of Saturn is dominated by the corotation torque, (c) Evolution 
of the inner/outer Lindblad torques and the corotation torque exerted on Jupiter. The 
average value of corotation torque vanishes and the migration of Jupiter is dominated by the 
Lindblad torques. The surface density profile is a = o~o r ~ 3 ^ 2 ■ (A color version of this figure 
is available in the online journal.) 
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Fig. 21. — This figure shows the cross section of surface density at different moments. 
The initial density profile is a = o"or~ 3 / 2 . The common gap forms after IOOOP10 evolution 
time. One may notice that the density gradient is positive ^ > at the out edge of the 
common gap, which produces positive corotation torque who drives Saturn outward away. 
The coordinates of planets in the figure denote their positions (x axis) and the average gas 
density around them (y axis). (A color version of this figure is available in the online journal.) 
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Fig. 22. — Panel (a) shows the evolution of mass variations within the coorbital zone of 
Jupiter and Saturn. Panel (b) zooms in from T = IOOOP10 to T = llOOPio- Panel (c) shows 
the evolution of corotation torque. Panel (d) zooms in from T = IOOOP10 to T = llOOPio 
too. The initial surface density profile is a = o§r~^l 2 . One may find the oscillations of 
the mass variation and the corotation torque match well. The short period of oscillations 
correspond to the orbit period and the long period is the libration period of the horse-shoe 
orbit, which equals to 32P 10 at the position of Saturn. (A color version of this figure is 
available in the online journal.) 




Fig. 23. — Density contours at T = 500Pio and T = 8OOP10 in a disk where a = <Tor~ 3 / 2 . 
One may find the Jupiter has opened a clear gap while the Saturn is still surrounded by gas 
at T — 500Pio (left figure). And the shock waves generated by Jupiter perturbs the coorbital 
zone of Saturn when the two gaps are overlapping T = 8OOP10 (right figure). 



